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A FINITE-VOLUME SCHEME FOR A SPINORIAL 
MATRIX DRIFT-DIFFUSION MODEL FOR SEMICONDUCTORS 

CLAIRE CHAINAIS-HILLAIRET, ANSGAR JUNGEL, AND POLINA SHPARTKO 


Abstract. An implicit Euler finite-volume scheme for a spinorial matrix drift-diffusion 
model for semiconductors is analyzed. The model consists of strongly coupled parabolic 
equations for the electron density matrix or, alternatively, of weakly coupled equations 
for the charge and spin-vector densities, coupled to the Poisson equation for the elec¬ 
tric potential. The equations are solved in a bounded domain with mixed Dirichlet- 
Neumann boundary conditions. The charge and spin-vector fluxes are approximated by 
a Scharfetter-Gummel discretization. The main features of the numerical scheme are the 
preservation of positivity and L°° bounds and the dissipation of the discrete free energy. 
The existence of a bounded discrete solution and the monotonicity of the discrete free 
energy are proved. For undoped semiconductor materials, the numerical scheme is uncon¬ 
ditionally stable. The fundamental ideas are reformulations using spin-up and spin-down 
densities and certain projections of the spin-vector density, free energy estimates, and a 
discrete Moser iteration. Furthermore, numerical simulations of a simple ferromagnetic- 
layer field-effect transistor in two space dimensions are presented. 


1. Introduction 

The exploitation of the electron spin in semiconductor devices is one of the promising 
trends for future electronics. Since the electron current can be controlled without changing 
the carrier concentration, this may allow for (almost) energy-conserving and fast-switching 
devices and, more generally, for electronic devices based on new operating principles. In 
the literature, several models have been proposed to describe the spin-polarized transport 
in semiconductor structures il |23]. Drift-diffusion approximations are widely employed 
[T8l |22] , since they do not require large computational resources but still describe the main 
transport phenomena. In this paper, we aim to analyze a hnite-volume scheme for a spin 
drift-diffusion system. Before we explain the model equations, we sketch the state of the 
art in spinorial drift-diffusion modeling. 

The existing drift-diffusion models can be classihed into two main groups. The hrst 
group is given by two-component drift-diffusion equations for the spin-up and spin-down 
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densities. One version of this model was rigorously derived from the spinor Boltzmann 
equation in the diffusion limit with strong spin-orbit coupling (compared to the mean-free 
path) [ 6 ]. A mathematical analysis of the limit model was performed in m, proving the 
global-in-time existence of weak solutions and their equilibration properties in two space 
dimensions. In three space dimensions, the well-posedness of the stationary system was 
shown in [10]. A quantum correction of Bohm potential type was derived in [3]. 

The second group consists of spin-vector drift-diffusion models in which the spin variable 
is a vector quantity. Combining the charge density with the spin-vector density, we can 
dehne the electron density matrix which solves a spinorial matrix drift-diffusion system. 
These models can be derived from the spinor Boltzmann equation by assuming a moderate 
spin-orbit coupling (B]. Projecting the spin-vector density in the direction of the precession 
vector, we recover the two-component drift-diffusion system as a special case. In | 6 |, the 
scattering rates are supposed to be scalar quantities. Assuming that the scattering rates 
are positive dehnite Hermitian matrices, a more general matrix drift-diffusion model was 
derived in [19]. The global existence of weak solutions to this model was shown in |16j . 

The aim of this paper is to analyze an implicit Euler hnite-volume approximation of the 
spinorial matrix drift-diffusion model of [TB] and to present some numerical simulations in 
two space dimensions. A numerical analysis of a hnite-volume scheme of the stationary 
two-component drift-diffusion equations was performed in HD]. A hnite-element scheme 
for a spin-vector equation with given electron current density (but coupled to the Landau- 
Lifshitz-Gilbert equation) was analyzed in [1] and simulated in [2] . However, no numerical 
analysis seems to be available so far for general spin-vector drift-diffusion models. 

1 . 1 . The model equations. The spin-vector model of [TB], which is analyzed in this 
paper, consists of the drift-diffusion equation for the (Hermitian) electron density matrix 
N G and the current density matrix J G 

1 

(1) dtN -I- div J + i'^[N^ m ■ a] = - \ -ii{N)ao — N 

T 

(2) J =+ NVV)p-^/‘^ inH, t>0, 

where [A, B\ = AB — BA is the commutator of two matrices A and B and H C is a 
bounded domain. The scaled physical parameters are the strength of the effective magnetic 
held, 7 > 0, the (normalized) direction of the precession vector m = { mi , m2 , m3) G the 
spin-hip relaxation time r > 0, and the dihusion coefficient H > 0. The precession vector 
plays the role of the local direction of the magnetization in the ferromagnet. In the analytic 
part of this paper, we assume that the precession vector m is constant. Furthermore, 
a = (cTi, ( 72 , CT 3 ) is the triple of the Pauli matrices (see [191 Formula (1)]), (Jq is the identity 
matrix in tr(A^) denotes the trace of the matrix N, and P = aQ + pm ■ a, where 

p G [0,1) represents the spin polarization of the scattering rates. The product m ■ a equals 
miCTi -|- 7712(72 + ^ 3 ( 73 . The electric potential V is self-consistently given by the Poisson 
equation 

(3) 



XjjAV = tr(A^) — C{x) in fl. 
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where Ad > 0 is the scaled Debye length and C(x) denotes the doping profile [H]. The 
bonndary and initial conditions are specihed below. 

In this paper, we investigate a scalar form of eqnations For this, we develop N 

and J in the Pauli basis via N = ^no<Jo + n- a and J = ^joao +j-a, where uq is the electron 
charge density and n the spin-vector density. Setting n = (rii, 77 , 2 , ns) and j = (ji,^ 2 ,^ 3 ) 
and defining r] = -^/l — , system 0-(§ can be written equivalently (see m Remark 1]) 

as 

(4) ^ino + divjo = 0, 

(5) dtUi + div ji - 2'j{n X m)i =- — , £=1,2,3, 

T 

D ^ D f ^ P \ 

jo = ^(-^0 - 2pJ ■ m), ji =— \^r]Ji + {l-r]){J-171)1711--Jomej , £=1,2,3, 

(7) Jo = —Vno — noVR, J = (Ji, J 2 , J 3 ) = —Vn — nVR in D, t > 0 . 

Moreover, the Poisson equation (|^ rewrites 

(8) — XjjAV = no — C{x) in D. 

System (|^-([^ is strongly coupled due to the cross-diffusion terms in ([^ and nonlinear due 
to the Poisson coupling. Note that any solution (no,n) to 0-0 defines a solution N to 
0-0 and vice versa. 

The boundary dQ = P'^ U P^ is assumed to consist of the union of contacts P^ and the 
isolating boundary part P^. Then the boundary and initial data are given by 


(9) 

no = n^. 

n = 11, V = 

on P^, t > 0, 

(10) 

Vno ■ V = Vni ■ v = W 

■ z/ = 0 on P^, t > 0, £ = 1,2,3, 

(11) 

’^o(O) = nl. 

n(0) = 

in D, 


where u is the exterior unit normal vector to dQ. 


1.2. Mathematical background. We aim to design a numerical scheme which preserves 
some qualitative properties of the continuous model, in particular preservation of the posi¬ 
tivity of the charge density, boundedness of the density matrix, and dissipation of the free 
energy. The main difficulty of the analysis is the strong coupling of the equations (the dif¬ 
fusion matrix is not diagonal), since maximum principle or regularity arguments generally 
do not apply. The key idea is to introduce two transformations of variables which make 
the diffusion matrix diagonal and thus reduce the level of coupling. 

The first transformation is defined by the spin-up and spin-down densities n± = |no ± 
n ■ m. Then system 0-0 becomes 

+ div (D(1 -|-p)(—V n+ — n+VR)) = —^(n+ — n_), 

dtn_ -I- div (iA(l — p){—Vn- — n_VR)) = — n+) 


( 12 ) 

(13) 
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and the boundary conditions 0. ( p!o| imply 
(14) 


n 


D 


n± = — on and Vn± ■ i/ = 0 on t > 0. 


N 


We observe that 0-0 implies ([T^-([T^ but not vice versa. Physically this is clear 
since the spin-up and spin-down densities contain less information than the full density 
matrix N. By the Stampacchia truncation method, the positivity and boundedness of 
n± was shown in [16], thus giving the positivity and boundedness of the charge density 
no = n+-|-n_. Using the notation a± = a+-|-a_, the (relative) free energy of the above 
system is given by the sum of the entropy and the electric energy. 


( 16 ) 


£( 4 ) = ^ J ^n±(\ogn± - 1) - n±\og^ + dx + ^ j \V{y - V'^)\^dx. 


Some formal computations show that it is nonnegative and nonincreasing for f > 0. 

The second transformation is given by the decomposition of n in the parallel and per¬ 
pendicular components with respect to m: n\\ = {n ■ m)m and n_i_ = n — n||. The equation 
for n_i_ reads as 


( 16 ) 


dtUy + div I — (—Vn_L — n^VU) ) — 27 (n_L x m) 
\V J 


T 


In [TB], it was proved by a Moser-type iteration technique that ny is bounded. Since 
n\\ = ^(n+ — n-)rh is bounded as well (see above), this implies an L°° bound for n and 
consequently for the density matrix N = ^nocxo + n ■ a. 

The task is to “translate” these ideas to a hnite-volume setting. We approximate the 
diffusive and convective part of the fluxes simultaneously by using a Scharfetter-Gummel 
discretization. These fluxes were introduced by Il’in [T3| and Scharfetter and Gummel 
[20] for the classical drift-diffusion model (without spin coupling). The discretizations are 
second-order accurate in space and preserve the steady states. The dissipativity with an 
implicit Euler discretization was shown in [9]. The discrete steady states were proved to 
be bounded [101 . Discrete entropy (free energy) estimates and/or the exponential decay 
of the free energy along trajectories towards the global equilibrium were investigated in 
[5l [I2] but still without any spin coupling. 

Our main results, detailed in Section are the existence of a bounded discrete solution 
to a fully discrete finite-volume scheme for ([4l)-([IT|) and the monotonicity of the discrete free 
energy for the spin-up and spin-down densities. The mathematical challenge is the proof of 
lower and upper bounds for the discrete densities. The “translation” of the Stampacchia 
truncation argument from the continuous to the discret case in, e.g., (16) faces some 
difficulties due to the drift term. The main difficulty lies in the fact that the monotonicity 
of the drift term (with respect to the density variable) cannot be exploited. Therefore, a 
Moser-type iteration method was employed in [16] . The idea is to derive a uniform estimate 
for n_[_ in the norm of the form 

< cg||n±||^. 


f > 0, 
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where || ■ ||g denotes the L'^{Q) norm and c > 0 does not depend on g e (l,cx)). By 
Gronwall’s lemma, this implies that 

\\n{t)\\g < e'=*||n^(0)||g, 

and the limit g —)■ cxd shows the claim. The discrete eqnivalent of the estimate is 

where n’]_ is an approximation of n± at time and At is the (nniform) time step size. In 
order to solve this recnrsion, we reqnire 1 — eg At > 0, thns imposing a condition on the 
time step size for fixed g. This motivates additional conditions on the model parameters, 
which are described and disenssed in Section 12.31 

The paper is organized as follows. In Section we detail the nnmerical scheme and 
present the main results, in particular the existence of discrete solutions (Theorem and 
the dissipativity of the discrete free energy (Theorem |^. The proofs are given in Sections 
1^ and 1^ Some numerical tests are presented in Section 

2 . Numerical method and main results 

In this section, we specify the numerical discretization of the spin drift-diffusion system 
0-0 and state the main results of the paper. 

2.1. Notations. Before we state the numerical scheme, we need to define the mesh of the 
domain hi and to introduce some notation. We consider the two-dimensional case only but 
the scheme can be generalized in a straightforward way to higher dimensions. 

Let n C be an open bounded polygonal set. The mesh A4 = {T,S,V) is given by 
a family T of open polygonal control volumes or cells, a family £ of edges, and a family 
V = {xk)k&t of points. We assume that the mesh is admissible in the sense of Definition 
9.1 in [ 7 ]. This definition implies that the straight line between two neighboring centers 
of cell {xk,xl) is orthogonal to the edge a = K\L between two control volumes K and L 
and therefore collinear to the unit normal vector to a outward to K. For instance, 
triangular meshes satisfy the admissibility condition if all angles of the triangles are smaller 
than 7r/2 [71 Example 9.1]. Voronoi meshes are also admissible meshes [71 Example 9.2]. 

Each edge a E £ is either an internal edge, a = K\L, or an exterior edge, a C dfl, 
and we set £ = Ant U Axt- We assume that each exterior edge is an element of either the 
Dirichlet or Neumann boundary such that we can set Axt = £^t For a given control 

volume K E T, we dehne the set £k of the edges of K, which can be written as the union 
of £K,mt, ^Kexti T^gxf For every a G T, there exists at least one cell K eT satisfying 
cr G £k, and we denote this cell by K„. When a is an interior cell with a = K\L, we have 
= K OT = L. 

For K E T and a E £k, we denote by dK,a the distance dK,a = d^xxyCr). Then, for 
a E Ant, cr = K\L, we define do- = dK,a + dL,o- = d^xxyXi) and for cr G Axt with cr G £k, 
do- = dx,o-- Furthermore, the measure of cr G T or a set a; C D is denoted by m(cr) or m(a;), 
respectively. In the numerical scheme, we need the so-called transmissibility coefficient 
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To- = ni((T)/do- for a E S. We assume that the mesh satishes the regularity constraint 

(17) 3^ > 0 : V77 G T : 'ia E £k '■ > ^diam(77). 

The hnite-volume scheme for a conservation law with unknown u provides a vector 
Ur = {uk)k£T of approximate values and the associated piecewise constant function, still 
denoted by ur, ur = which approximates the unknown u. Here, Ik denotes 

the characteristic function of the cell K. The approximate values of the Dirichlet boundary 
provide a vector u^d = The vector containing the approximate values in the 

control volumes and at the Dirichlet boundary edges is denoted by um = {.ur,U£D). 

The numerical scheme can be formulated in a compact form by introducing the following 
notation. For any vector Um = {ur,U£D), we dehne, for all iF G T and a E £k, 


Ur, a — 



if a = K\L, 

if 

•r _ nN 

li a — 


and we set T)uk,u = — ur- We remark that the dehnition of ur^^ ensures that 

D'Wx,o- = 0 on the Neumann boundary edges. Then the discrete seminorm for um can 
be dehned by 

1/2 


\um\im — I 

VcrSf 


d|Dm 


K,(t\ 


where the summation is over all edges a E 8 with K = K„. The U’ norm of ur reads as 

i/p 


\M\p = 1 m(iF)|Mi^|^’ 

yKeT 


for 1 < p < oo and ||nr||oo = max \ uk \. 


When formulating a hnite-volume scheme, we have to dehne some numerical huxes Jx,cr 
which are consistent approximations of the exact huxes through the edges J ■ UR ^jCts. 
We impose the conservation of the numerical huxes Jx,(t + JL,a for = K\L, requiring that 
they vanish on the Neumann boundary edges, Jr^o- = 0 for ct G Tj^ext- Then the discrete 
integration-by-parts formula becomes 


(18) 


EE JK,aUK — 

K&T o-GSk 




ct££ 




2.2. Numerical scheme. At each time step k > 0, we dehne the approximate solution 
Ur = {u^R)K^r for u G {no,n, H} and the approximate values at the Dirichlet boundary. 


- K'^aJa&S, 


= [U 


(u^)^^£d (which in fact does not depend on k since the boundary data is time- 


independent). We hrst dehne the initial and boundary conditions corresponding to 0 
and (^. We set 


n, 


0 


?o 


0,A) "•A 


n'l' = 


m(iF) 


(ng,n“)(ix for all K eT, 


>K 


(19) 
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= —TT (n , 0 , 1 / )ds for all a G 
m(a) 

Note that = 0 for cr G We may define similarly the quantities Ck, Dk,Pk for 
a given K E T. 

We consider a temporal implicit Euler and spatial hnite-volume discretization. The 
scheme for Q, ([^, (|^ writes, for all it' G T and A; > 1, as 


( 20 ) 

( 21 ) 

( 22 ) 


m(it') 


m(it') 




At 


jo,K,a — 0, 

(t&Sk 


TI'^-TIk 


k-1 


At 


+ Y1 - 27 m(it:)(n^ x = 

(tGSk 


m(K) 


-n 


T 


A, 


-XlY. = m{K)(nl„ - Ck), 

o-&£k 

where the discrete counterpart to (| 6 |) is, for all it' G T, a G Ek, k >0, 


(23) 

(24) 


jo,K,r7 — 2 , 

1(7 




dK,a = -f + (1 “ V^)iJK,a ' - '^Jo,K,a^c 

Va ^ ^ 


7k 


P<r jk 


The numerical fluxes approximations of the integrals Jo ■ ^K,ads 

and Ji ■ i'K,ads at time kAt, and we set = (J^x,(t)^=i, 2 , 3 - recall that Jq and J 
are dehned by Q. We use a Scharfetter-Gummel approximation for the dehnition of the 
numerical fluxes. For given it' G T and a G Ek, we set 

(26) 4K,, = T4B{DVl,)nlK-B{-DVl,)nl.KP, « = 0,1,2,3, 

where B is the Bernoulli function dehned by 

X 


B{x) = 


exp(x)— 1 


for X 7 ^ 0 and B{0) = 1. 


It remains to dehne the quantities rho-, Pa and appearing in (23) and (24). We 
use a weighted harmonic average on the interior edges and a classical mean value on the 
boundary edges. 


= 


d(7 DkDl 


dA,cr Dl + dLL,a Dk 


for a G £^int, cr = K\L, D„ = 


m (T 


D{s)ds for a G 


and similar dehnitions for ruo- and p^- Furthermore, we set — p^. 

Finally, the boundary conditions are 

(26) = n^ = 0, = for a G 

(27) UnlK., = = 0 for e C = 0.1, 2, 3, i > 0. 

We remark that they imply = 0 for a G ^ = 0,1, 2, 3, and fc > 0. 
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For later use, we note that, using the elementary property B{x) — B[—x) = —x for 
X G M, the numerical fluxes can be reformulated in two different manners: 

(28) JIk., = r,( - - B(-DVt,)Dnp.>) 

(29) <’ = 0.1,2,3, 

and adding these expressions leads to a third formulation: 


(30) 
where 

(31) 


— 'o- 


-2Ka + - B\I^Vl^ 


)Dn: 




B^ 


, , X , /x\ Bix) 

(U=2C0th(-)=-U 


B{-x) 


2.3. Main results. We impose the following assumptions on the domain and the data: 


(32) 

(33) 

(34) 

(35) 


C bounded domain, dVt = F'^ U F^, 
D, p, m are constant and |m| = 1, 


pD n F^ = 0, m(F") > 0, F^ open. 


no, n 


G 


-Uq ± n° ■ m > 0, > 0, n^, G 


M = {T^£,V) is an admissible mesh satisfying (17). 


We hrst remark that if n^, V^) is a solution to scheme (20)-(27) for a given k>l 
((uQ-y-, n!^) are dehned as the discretization of the initial conditions), we can dehne n^^ = 


^Uq j- ± ■ m, n\-Y = fij — (n^ ■ m)m. Moreover, as and are dehned on the whole 

domain f2, we can dehne and V-y by taking the mean value of and on each 
control volume K E T. 

Then the following existence result holds. 

Theorem 1 (Existence of a solution to the numerical scheme and L°° bounds). Let as¬ 
sumptions (p^-(35) hold. We impose the following constraints: 

pXj 


(36) 


1 \2 

0 < At < - ■=__ 

-a' C(l+p)||C|L 


0 < r < 




mcl 


Then for k > 1, there exists a solution to scheme (20)-(27) satisyfing 


0 < < 2M°, 0 < < M' 


n. 


ri - 


< 2M^ in n, 


where = M°(l — aAt) ^ and 


= max ( - sup n 
2 an 


D 


sup 

o 


-Uq + \n ■ m\ 


sup |u 


0 I 


sup C 
n 


In the continuous case, similar L°° bounds for the spin-up and spin-down densities, and 
therefore for the electron charge density, were shown in [16]. These bounds do not depend 
on time. The mixing of the spin-vector components prevents the use of the monotonicity 
argument for nj_, solving (16). Therefore, both in the continuous and discrete situations, 
the bound for the spin-vector density depends on time. 
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The constraint on At is needed in the definition of M^. Fnrthermore, the condition on 
r is necessary to prove the L°° bonnd for n\ We believe that the latter restriction is 
technical. We stress the fact that onr scheme is nnconditionally stable if the semicondnctor 
is nndoped, i.e. (7 = 0. In this situation, At and r can be chosen arbitrarily. 

Next, we prove that the scheme dissipates the discrete free energy, defined by 

= E E - 1) - <K log ^ + ’f) 

(37) -'''“)*>)"■ 

(j^S 


Theorem 2 (Dissipation of the discrete free energy). Let assumptions (32)-(35) hold and 
let (uq . 7 -, V^)fc>o be a solution to scheme (|20|)-(27) satisyfing 0 < 11 ^ 7 - < M°. We 
further assume that > n* > 0 and that log(n^/2) -|- is constant in D. Then the 
mapping k ^ 


At 


(38) B‘+ — 


is nonincreasing, i.e., the scheme dissipates the free energy (|37|).' 

{4./t."W}(D(log4 + U)K,,)" < E 


± 


(T&E 


T„ mini 


nfc—1 


k>l. 


The above dissipation inequality for the free energy is the discrete counterpart of the 
continuous estimate for the free energy (15) [TBl Formula (28)]: 


E{t) + - 



0 Jn 




D{1 ±p)n±|V(logn± + V)\‘^dxds < E{0). 


One may ask if the discrete solution converges to the continuous one when the approx¬ 
imation parameters tend to zero. However, it seems to be difficult to extract a discrete 
gradient estimate for 7 - from the discrete free energy estimate in Theorem since we 
do not have a suitable discrete version of the chain rule n±|Vlogn±p = 4|Vyd7±P- 


3. Proof of Theorem [T] 

The proof of Theorem [T] will be presented in two subsections. We first establish the 
existence of a solution (uqt-, F^) at each time step fc > 1 by an induction argument. 
The proof is based on the fixed-point theorem of Brouwer. In this subsection, we also show 
L°° bounds on Uqj- and 7- which depend on k. Then, in the second subsection, we prove 
that these bounds are in fact uniform with respect to k. 


3.1. Existence of a solution to the scheme. We first note that the initial condition 
(rig 7 -, n^) is well-defined by (19). Moreover, the definition of ensures that In’} 7 -I < M°, 
0 < 7 - < and therefore 0 < rinT- < 2M° and \n^ -t\ < 2M°. 


^ "•o,r — 




^k—1 zfk—1 


The proof is done by induction. Let k > 1. Assuming that 


,VJf ^) is given 


and verifies ^ 


< ^ we will prove the existence of (uqt-, n^, VJf), 

solution to (20)-(27), satisfying these bounds with k instead of k — 1. Scheme (20)-(27) is a 
nonlinear system of equations. We prove the existence of a solution by using a fixed-point 












10 


C. CHAINAIS-HILLAIRET, A. JUNGEL, AND P. SHPARTKO 


theorem. Let us denote by 6 the cardinality of the mesh T (the number of control volumes) 


and let > 0. We dehne an application 


1)40 


-)■ 


1)40 


such that (pr) = where 


pr = (po^t^Pt) and nr = (no,r)'^r)- If is based on a linearization of the scheme and 
dehned in two steps: 

• First, we dehne Fr G as the solution to the linear system 

Y,r^^VK,. = m{K){po,K-CK) for KeT, 

uGSk 


(39) 


■^D 


(40) 

(41) 


14 = for a e BVK,a = 0 for a G 

* Second, we construct nr = (no,r, ^t) ^ as the solution to 

miK) . miK) . . v-^ . , ^ 

(Ro,a — 'i1o,a) + P (^o,A — Po,k) + 2_^ jo,K,cT = 0 for iF G T, 

(t&Ek 

m(iF), m(iF) , , . 

) + p [ni^K - PpKj + 2^ Je,K,a 

o-GSk 


\/^ m(iF) 

27 m(A)(nK x m)e = - ne^K 


for iF G T, i = 1, 2, 3, 


where jo,A,cr and je,K,a are dehned in (23) and (24), with Je^K,a dehned in (25), but 
without the superindex k. The boundary conditions read as 


no, a = nf , n^ 


= 0 for (T G 4 


D 

ext’ 


= 1,2,3, 


= 0 foraG^^.,,,, iFGT, £ = 0,1,2,3. 


■'A',ext 


(42) 

(43) Bn^^K.a 

The parameter p > 0 allows us to prove unconditional stability for the linearized prob¬ 
lem; see e.g. [1]. The corresponding term vanishes for hxed points pr = nr, so that a hxed 
point for is a solution to scheme (20)-(27). We choose 


(44) 


P > 


mcl 

>^1 




1 1+p 


At. 


The existence and uniqueness of Vr, solution to (39), are obvious since the corresponding 
matrix is positive dehnite. As this matrix does not depend on pr and the right-hand side 
is continuous with respect to pr, the hrst mapping P 7 - i—>■ VV is continuous from to 
M®. This property is not so obvious for the second mapping, based on the linear system 
of equations (40)-([4^. We will prove this property below (Step 1), in order to guarantee 
that the mapping is well-dehned and continuous. 

Then, in order to apply Brouwer’s hxed-point theorem, we will prove that Fj^ preserves 
the set 

(45) 


— {r-t — iFopr-iTir) G : 0 < 4 Ir-UjtI 4 . 


It is a bounded set because each element nr G verihes 0 < no,r 4 2M^ and \nr\ 4 2M^. 
This part of the proof is the most challenging one. Given pr G and nr = F4pr)) we 
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will first establish the positivity of n±j- (Step 2 ), then the upper bounds for n±^j- (Step 
3), and hnally the bound for n±^j- (Step 4). 


Step 1: Existence and uniqueness of solutions to (40)-(43). The linear system of equa¬ 
tions (40)-(43) is a square system of size 46. The existence of a solution is equivalent to the 
uniqueness of a solution and to the invertibility of the corresponding matrix. Therefore, 
we just have to prove that if the right-hand side to the system is zero then the solution is 
zero. Thus, we may work with the original linear system assuming homogeneous Dirichlet 
boundary conditions and setting = po.A = 0 and = pW = 0 , in order to set the 
right-hand side to zero. 

We multiply the corresponding equation (40) by and (41) by ne^x, sum these four 

equations, and sum over all control volumes K E T: 


0 - (1 + p) ^ x + {1 + ^ + 4 X] X] jo,K,ano, 

KeT KeT KeTaeSK 


K 


K&T o-eEji 
= Ti + ■ ■ ■ + Tg. 


27 E m{K){nK X m) 


m { K ) ^ 

UK + 


AST 


r 




AST 


Note that T 5 = 0 and Ti, T 2 , and Tg are nonnegative. Thus, it remains to estimate the 
terms and T4. By discrete integration by parts (note that the problem is homogeneous) 
and the dehnitions ([^-(24) of j(,^K,a (omitting the superindex fc), 

T-i = —7-7J ^^(^o,A,o- — ‘^pJK,a ■ 'ih)Dno,A,o- =: T31 -|- T32, 


T - ^ 

J-i — — 5 - 

jjZ 


tGS 


(vJK,a + (1 - v){JK,a ' m)m - ' DnA,<7 =: Tn -F T42 -h T43. 

(tGS 


With formulation (30), dehnition (31) of B^, and the discrete chain rule 
= D(n^)A,(T, we have 

T 31 = { 2 B\J 4 VK,a){E)no,K,af + E){nl)K,aE)VK,a) , 


7^8 


T32 — 


pD 

At]"^ 


r^(2S*(D14',<T)DnA,f7 ■ m + [uk + nK,a) ■ ThDVK,a)Eino^K,a, 


tG8 


D 


^41 = — r,, {2B^{DVK,a)\BnK,a\‘' + , 

T 42 = 5^ry25^(DVK,,)(Dr?,,,^ ■ rhf + D{{n ■ mf)K,aBVK,a), 


T43 = 


2^2 

pD 


T^S 


477 ^ 




a£S 
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We collect all terms from T 3 + T 4 involving the fnnction 




TaB'^{DVK,a)(^(Dno^K,aY + ri\DnK,a\^ + (1 - v)(PnK,^ ■ mf 

ae£ ^ 


pDno^K,a'DnK,a ■ rn 


D 




( Di 1 o,A,o 


7^S 


1/4 -p/2 


I Dni^^o- ■ hi, j V —p/2 1 — 1/^/2 j \ Dfix,cr ■ hi 


Drii 


0,A,cr 


2 

+ y ^ ~ 2 ) • rri|^) 

The eigenvalnes of the 2x2 matrix appearing in Ji are 

A± = ^(5 - 2r/2) ± iV(5-2r/2)2-8r/2 > ^ > 0. 

Then, nsing the ineqnalities B^{z) > 1 for all x G M and \Y)nK,a\^ > \DnK^a - hip (since 
|mp = 1 ), it follows that 


-fi > |^A_((Dno,i^,^)^ + {DnK,a ■ mf) + y|Dhx,<^|^ ) > 0. 


Next, we collect in I 2 the remaining terms from T 3 + T 4 involving the discrete gradient 
DVr-.ct- Taking into acconnt that 

{uk + nK,a) ■ fnDnQ^K,a + {no^K + ?^o,A,f7)Dh/f,<7 ■ hi = 2D((h ■ hi)no)A-,a, 


integrating by parts, and employing the discrete Poisson eqnation (|39|), we infer that 
D 


h : = 


2 p 


f + viY){\n\^)K,a + (1 - p)D((h ■ hi)^)^,, 


pD((h • m)no)K, 
D 

(T 

K&T 


2r]^Xl 


T^i^)iPo,K - Ck) f + (1 - ■ hi)^ 


- p{nK ■ rn)no,K 


The snm of the terms in the brackets is nonnegative since 


+ vI^kI^ + (1 - v){nK ■ mf - p{nK ■ m)no,K 
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Therefore, the term involving po,A ^ 0 can be omitted, giving 
D 


h> 


> 




D 


2r/2A^ 


K&T 

This shows finally that 


m{K)CK f k + (1 - - Pi^K ■ 'rh)no^Kj 

K£T ^ ' 

IlC'lloo m(ir) {^nl^ + = 




liq|oo(^||no,r||^ + ||nr||^ 


T4> 


mcl 


D 


-J^noAl + WnrWl^. 


and summarizing the estimates for Ti,..., Tg, we conclude that 

T + /i D\\Cl 


At 




^ll’^o,rll2 + Il^rll2^ < 0. 


Hence, choosing /i as in (44), the first bracket is positive, showing that no,r = 0 and fij- = 0, 
which proves the invertibility of the linear system of equations (40)-(43). The second step 
involved in the definition of (VV, Pr) n-j-, is a well-defined mapping. Moreover, 
the matrix and the right-hand side of the linear system of equations are continuous with 
respect to (VV, pr) so that the mapping is continuous. 

Step 2: Positivity of n±,r- We will prove that n±^K > 0 for all K E T. Multiplying 
(40)-(41) by m and adding or subtracting it from (40), multiplied by we hud that 


(46) 


- 4^) + + D(l±p) ^ 


At 


At 

m{K) 

= — {n+,K-n_^K), 

It 


o-gSk 


where p±^K = \po,k ± Pk ■ rfiK and J±,K,a = \Jo,K,a ± Jr, a ■ m, i-e. 
(47) J±,K,a = Tu{B{Y)VK,a)'n±^K “ B{ — 'DVK,a)n±^K,a)- 


Then, multiplying (46) by ^ = min{0, summing over all control volumes K E T, 

and adding both equations, it follows that 


\^ \^ m(77) , k—^ \ — \ —^ \^ m(77) , . _ 

0 = -^{n±,K - “ P±,k)^±,k 

± AST ± AST 


(48) 


+ -O 5^(1 ± p) ^ J±.K.,n^j! + A 

± K£Ta£SK KeT 

=: Tj + Tg + Tg + Tio, 


m(F) 
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Since ^ > 0, and p±^K > 0, the first two terms in (48) can be 

estimated as 

± AST ± Aer 

The monotonicity of the mapping z ^ z~ shows that Tiq is nonnegative. By the discrete 
integration-by-parts formula (|T^ , the third term in (48) becomes 


Tg = -D ^(1 ± p) ^ J±,K,a'D{n^)K,o 




The sum over the boundary edges vanishes since = 0 for all a G We claim that 

(49) - J±,A,aD(n^)i^,^ > - {u^^kT), forKeT and a e Ek, 


such that 

( 50 ) 


D 


n>-'£(l±p)'£ T„DV^,,D((nif) 


K,a' 


ct££ 


To prove (49), we distinguish the cases DVk,(t > 0 and DVk,(t < 0. If T^Vk,^ > 0, we 
apply a formulation similar to ( [^ , leading to 

'ht,A,(TD(7i^) a,(t Per(Df^jO-R-itiAjcr T B(JDVj^u^Tyn^ j^u^E)Xij. 

Then, using the nonnegativity of the function B, the monotonicity of the mapping z ^ z~, 
and the inequality 2 ;( 2 ;“—i/“) > — {y~)'^), we obtain (49). If DV/^o- < 0, we employ 

formulation (28), so that 

~'^±,A,crD(7i±)A,cr = Eo-(Df^A,crR-±,A + -B(—Df^A,(T)D7i±^x,cr)D(?^,^)x,<T 
and similar arguments lead to (49). 

Applying discrete integration by parts to the right-hand side of ( |50| ) (the boundary 
term vanishes since the boundary data is nonnegative) and employing the discrete Poisson 
equation (39), we hnd that 

Tg > -y ^ 7-<7DlAf,<,((l+p)(n"^)^ + (1 -p)(n:^^)^) 

AeT o-gSk 

= 7 ^ 5^m(A:)(po,A-C';^)((l+p)(n“^)2 + (l-p)(n:_^)2) 

Ker 


> 


D 


2X1 


IlC'lloo +P)(77 +,a)^ + (1 -p)(77_,x)^)- 


AeT 


Summarizing the above estimates, we conclude from (48) that 

(^ ■ E m(A-)((n;,^)= + («:,,)’) < 0. 

V D / 
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By the choice of jj, in ( [44| ), we deduce that 

+ (n-^)2) < 0, 

KeT 


which implies that = 0 and hence n±,_fc > 0 for all K E T. 


Step 3: Upper bounds for n±,r ■ The goal is to show that n±^K < for all K E T, where 
is dehned in Theorem]^ We multiply (46) by {n±^K — = max{0,n±,x — 

sum over all K E T, and add both equations: 


± AST 

+ f-Z E - f''^'‘))E±,K - MY 

± KeT 

(51) + E E - MY 

± KeT 

+ DYi^±P)J2Yl J±,KAn±,K - 

± K&Tu&Sk 

+ iEE m(i^)((n+,^ - M^) - {n_^K - M^)){±{n±^K - 

± KYT 

=: Til + Ti2 + Ti3 + Ti4 + Ti5. 


Using the inequality (z — y)z~^ > A(z~^)^ — the hrst two terms are estimated by 

L. > ^E E n'(A')((n±,K - MYf, 

± KeT 

^■2 > ^ E E YK){{n±.K - mY)\ 

± KeT 

since n^K ^ and p±^K E by assumption. By dehnition of M^, the third term 

Ti 3 becomes 

Ti3 = EE m{K){nYK-MY, 

± AST 

and the last term T 15 is nonnegative. 

It remains to estimate T 14 . By discrete integration by parts (the boundary term vanishes 
in view of for a E we find that 

Tu = -D 5^(1 ± p) A,K,,D((n± - MY)„,,- 

di a^£ 
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Similarly as in Step 2, we claim that the following estimate holds: 

ri4 > D 5^(1 ±p)Y, - M")+)2 + M\n^ - M^) 




K,a 


Indeed, let hrst T)VK,a > 0. Using the ineqnalities D(n± — — M^)^)K^a > 0 

and {n±^K,r 7 — M^)D((n± — M^)^)K,a > |D(((n± — M^)^Y)K,a > 0, it follows from (29) 
that 

-J±,K,. > r.DUx,.((n±,^,. - M^) + M'=)D((n± - 

> r.DU^,, QD(((n± - + M"D((n± - . 

The proof for DVx.a < 0 is similar, employing formnlation (28). Then, integrating by parts 
and employing the Poisson eqnation and po,A > 0, 

Tu > § E6 ±P) E “(/0(Po,*- - Ck) (^(("±,a- - M*)+)= + M‘(n±,K - M‘)+j 

> -|;l|C||o„(l + P) E E - J\A)+) . 


± K&T 

Snmmarizing the above estimates, we infer from ( [5T| that 
1 + /i D 

‘ ± KeT 


2At XI 


- ^iiciu(i+p)) E E “(^o((«±,A - A/‘)+)' 


a 


D 


liqu(i + p))M^ EE m(K)(n±,j^-M'^)+ < 0 . 


± AST 


Then, choosing jj, as in (44) and taking into acconnt the dehnition of a, we infer that 
ii±,A < for K E T. 

Step 4- L°° bound for n\ q-. We prove a nniform bonnd for = hk — {uk ■ rh)m. 
For this, we multiply the vector version of (41) (omitting the superindex k) by rh twice, 
and taking the difference of the equations for nx and {riK ■ m)rh, we obtain 


- n^±^K + Kn±,K - P±,k)) + — ^ J±,K,a - 27m(iF)(n_L,ii- x m) 

^ <t&£k 


At 


m{K) 


T 




where p±^K = Pk — {pk ■ m)m, and J±,K,a is given by 

(52) J±,K,a = a{ — DVK,crn±,A — B(—DVK,a)^n±^K,a) 

(53) = Tcr{ — ^VK,an±^K,a — B (DVK,a)^n±^K,a) ■ 
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Then, multiplying this equation by (where q E N) and summing over 

K E T, we arrive at Tig + Tn + Tig + Tig = 0 with 

Tie = ^ ^ m{K){n±^K - ^^ 7 ^) ■ n±^K\n±,K\‘^^'^~^\ 


At 


AST 


Ti 7 = ^ ^ m{K){n±^K - P±,k) • n±^K\n±,K\‘^^'‘ 


At 

D 


K&T 


TlS — — J±,K,a • n±^K\n±^K\‘^^'^ 


V 


K&T o-^Ek 


Ti9 = - ^ m{K)\n±^K 


\2q 


KeT 


The elementary inequality ■ (a — b) > (|ap'^ — |&p'^)/(2g) for a, b eM.^ shows that 

1 


Tie > 


Ti7 > 


2qAt 

P 

2qAt 


m(T:)|nl,K 


-1 \2q 

k\ 


\K&T 


K&T 


m(T')|n_L,K|^'^ - m{K)\p±^K\^‘^ 

K&T 


\K&T 


By discrete integration by parts (observe that f?_L^o- = 0 for a G £^t), 

Tig =-^,A,o- ■ {n±,K,a\n±,K,a\‘^^'^~^'' —'n±,K\'n±,K\‘^^'^~^^) ■ 

^ a&£ 

Again, we distinguish the cases DV/^ > 0 and T2)Vk,(t < 0 for given K E T and a E £k- 
First, let DV/^ cr > 0 and use formulation (53) of the numerical flux. This gives 

— ■J±,K,a ■ {n±,K,a\n±,K,a-\‘^^'^~^^ — n±^K\n±,K\‘^^'^~^'') 

= TrjDVK,a'n±,K,a ' {n±,K,a\n±,K,a\'^^‘^~^'^ — T_L^ii' 

+ T(DV7c,cr)(T_L,ii-,o- — ut.a) • {n±,K,a\n±^K,tT\‘^‘''^~^^ — 

Because of 

a- (a|a|2(''-') - = |ap'? - a ■ - (l - \b\^'^ 

2q \ Zq J 

> (l - ^) (lap-? - |6|29) for all a, 6 G 

applied to a = n±^K,a and b = n±^K, and the monotonicity of the mapping a h-)■ 
we hnd that 

— J±,K,(t ■ {n±,K,a\n±^K,a\‘^^'^~^^ — |n_L^x 

> To- D1 /x,o-(|Tt,A,oP'^ — \n±,K\^'^)- 
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This result still holds if < 0, thanks to formulation (52). Therefore 


D 


V 


Ti 8> — [I - — j y^^TaDVK,aB{\n±\ 


2 g' 


K,cr- 


t&£ 


Using discrete integration by parts and the Poisson equation ( |^ leads to 

D 

2q 


Ti8 > 


riXl 


D 


1 - Ck)W,k\^ > - JrllCIloo m(A')|.u,Af 


Aer 




Aer 


Summarizing the above estimates, we obtain 

D\\C\V 


+ /i + 


rjX 


D 


^ m(iP)|n_L,Kp'' < ^ m{K)\ 


n 


k-l \ 2q 

±,k\ 


K&T 


K&T 


fi m{K)\p±^K 


\2q 


K&T 


Condition (36) on r, the induction hypothesis ||u 5 ^ 7 ^||oo < ^ and the fact that 


p G (see (45) for the definition of 5^), such that ||pj 


,rl 


< imply that 


||if±,r|| 2 q < meas(r 2 )^'^*'^^^M^ for g > 1 . 

Passing to the limit q —>■ +cxd, we deduce that ||n_L,r||oo < 

Conclusion. In Step 1, we have proved that the mapping is well-defined and con¬ 
tinuous. In Steps 2-4, we have proved that preserves the bounded set S^. Thus, the 
hxed-point theorem of Brouwer shows the existence of a fixed point to F^, belonging to 
Let us denote this hxed point by = (uqt-, u^). It is a solution to scheme (20)-(27) 
at step k and satishes 

and |n^ ^| < for iC G T. 


,AI ^ 


3.2. Uniform bounds for the spin-up and spin-down densities. In order to conclude 
the proof of Theorem it remains to prove that the upper bounds on the spin-up and 
spin-down densities in fact do not depend on k. The positivity of these densities is already 
proved above. 

We assume as induction hypothesis that for all iP G T (this property is 

ensured for fc = 1 by the definition of M°). Scheme (20)-(21) implies that 


(54) 


m(iP) , 


n 


,fc- 

±,A 


J) + r>(i±p) = t 

o-gSk 


m(iP) 

2 r 


[n. 


,A 


n 


-,K)- 


As in Step 3 above, we multiply (54) by (n^ ^ — M°)+, sum over all iP G T and add both 
equations. This yields S'! -|- S '2 -|- S '3 = 0, where 

= U U - V”) - - my, 


± K&T 
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Sj = D 5^(1 ± p) JIkA<k - M'Y, 

± KeTaGSK 

Sa = ^ E - <k) (i<.K - - M")+). 

^ KeT 

It is clear that 5*3 > 0 and, by the induction hypothesis, that 

± KgT 

The term S 2 is the analogue of T 14 . Following the same ideas as in Step 3, we obtain 

S 2 > ^ E(1 ± P) E - Ca) (h(4,A - - M»)+) . 

But, as the positivity of and and the dehnition of 

ensure that no,A ~ ~ — Ck > ~ leading to S 2 > 0. 

Therefore, we infer that 

EE“W («a-m“)+)To, 

± KeT 

which yields the expected result. 


4. Proof of Theorem [2] 


Let V^)k>o be a solution to (22), (54) with the corresponding Dirichlet-Neumann 

boundary conditions. Since we have to deal with the logarithm of the densities which 
may vanish, we introduce a regularization of the discrete free energy. For 5 > 0, we set 


k,5 

^±,K = ^±,K 


+ 6 and dehne 


(55) yy 


± K&T 




k,5 

±,K 


) - 1 ) - ni^K log 




v^)K,ay. 


tGS 




n 




Therefore, we have ^ = Ui + U 2 , where 


Ui = '^Yl ^iK)(nX%{lognX% - 1 ) - nX}’\lognX}'^ - 1 ) 

± AST ^ 

- (ni^K - nX~K) log + '^) ) ’ 

U2 = ^ - V^)x,af - - V^)x,ay). 

ctGE 
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The convexity of a; h->• a;(loga; — 1) shows that a;(loga; — 1) — y{\ogy — l)<{x — y) logx for 
all X, ?/ > 0. Hence, 




logn^’Jf^ - log ( ^ + <5 


n 


'^1 s E E 

± AST 

Using the elementary inequality — y"^) < (x — y)x for all x, ?/ G M, integrating by 
parts, and employing the discrete Poisson equation (|2^, it follows that 


U2<XlYl - V^)K,a 

= E E - n? 


KeT o-gSk 


^ ^ - niy)(ri - v°). 


± AST 


We summarize the above inequalities and use scheme (|54|) to hnd that 
1 

Ai' 


- El ~ ^-,k) (- log’^-i) 

^ AST 


5^ D(i ± p) 5^ (log n‘y + rj - log (qt+i 

± K&T(7&£k 


D 


V, 


The first term on the right-hand side is clearly nonpositive. We apply the discrete integra- 


tion-by-parts formula (18) to the second term. Then, with the hypothesis on the boundary 
data (i.e. log(n^/2) -|- is constant in H such that = —D(logn'^)x,o- for all iP G T 

and cr G £^ic), we infer that 

i(i5f - Bp') < 5^ £>(1 ± p) 5^ J‘,,,,„D(logny + V'‘)k,, 

di 

+ 4,A,aD( log - log(n^ + 26))^^^. 

di (tG5 

Introducing the numerical fluxes associated to the regularized densities, 

we can write 

^(Bf - Bp') < Bs + f/4 + f/, + Be, 

= D{1 ±P)Y, Ti,,D( log"i' + 

it (tG5 

U^ = Y^ D{1 ±p)Y^ 4;'5^ ,^D(logn^ - log(n^ + 26))k, a, 

di (tG5 


where 
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± crg£ 

Ue = 6j2D{l±p)J2 logn'' - log(n'' + 25))^^^ 


di cr£5 

Now, we employ the following inequalities, which are proved in [U Appendix A]: 

'^±’A,aD(logn±’^ + V^)K,a < 

< r^max(n^’|^^,n^’Jj^_^)|D(logn±’^ + V^)K,a\- 

The first inequality yields 

f/3 < -^T>(l±p)^r^min(n^’|^^,nJ^’|^;^^^)(D(logn±’‘^ + l/^)i^,^)^ 
di a^E 

while the second one, together with Young’s inequality, gives t/4 < [/41 + t/42, where 

1 X ^ V / - \ X ^ . / h ^ h ^ \ / /■. t’A -r T-U\ \2 


f/41 = ^ Y min(n^y, (D(logn±’'^ + V^)K,a) , 

di (tG£’ 

Tj dm 4- 1 x\2(D(logn-^ — log(n-° + 25))ii'o-)^ 

t^42 = 2^T>(l±p) 2^r,(max(n±_^,ni’^^^)) - ^ -’ 

ib (tG5 


±,A) '^±,K,a) 


^ I log(n;;\^+ 5)-logn;^|^_^, 

since min(n^’'^;^, > h for all A' G T and a G £k- Applying Young’s inequality again, 

we obtain < t/51 + U 52 with 

U51 = t ^T)(l ±p) ^r^min(n^y,n^y^)(D(logn±’^ + 


(tGE 






(DVJ. 


and 


• / A:,5 /c,5 N 

ef m>>nn±>.n±>,,) 


Ue < 2Z)U I log(n° + 26) - lognj, 


<2Z)i5^r4Dl/,t 


ct££ 


Summarizing the above inequalities, we deduce that 
’6) -^{Es - E^-^) + ^ ^T)(l ± ni V T^mininb'^ 


7^S 


< AD5 Y,rAWl^)^ + 2D[6 + 


cr&£ 


(MO + 5)^ 


a£E 


k^S 

^±,K,c 

J(D(log^ 

log(n^ + 25) 


.D |2 


M' 


On the one hand, the term Xlo-ef does not depend on 5 and is bounded (th 


IS 


can be seen by using scheme (22) and the L°° bound on n^j-). On the other hand, we 
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rewrite 


log(n^ + 26)- logn, 


D 

M\ 


|2 

. ( 2(5 \ 

\iM ~ 



l,M 



log 



2 


Employing the inequality | log?/ — loga;| < \y — x\l min(a;, y) for x, ?/ > 0, and the fact that 
> n* > 0, we obtain 


\og{n^ + 26) - \ogn^\\ j^ < ^ 


n 


D |2 


Thanks to hypothesis (33), G H^{VL)^ and Lemma 9.4 in [7], we conclude that 
7L||n^||j:/i(Q) with K depending only on the regularity of the mesh Ai. Therefore, the right- 




hand side in (56) tends to zero when 5 —)■ 0. Passing to the limit 5 —>■ 0 in (56) then leads 
to ([M|). This concludes the proof of Theorem]^ 


5. Numerical simulations 

As an illustration of the numerical scheme, analyzed in the previous sections, we present 
two-dimensional simulations of a simple double-gate ferromagnetic MESFET (metal semi¬ 
conductor field-effect transistor). This device is composed of a semiconductor region which 
is sandwiched between two ferromagnetic contact regions (see Figure [^. The idea of 
such devices is that the source region plays the role of a spin polarizer. The non-zero 
spin-orbit interaction causes the electrons to process during the propagation through the 
middle channel region. At the drain contact, only those electrons with spin aligned to the 
drain magnetization can leave the channel and contribute to the current flow. Here, we 
focus on the feasibility of our numerical scheme and the verification of the properties of 
the numerical solution and less on the physical properties. Therefore, the physical setting 
considered here is strongly simplihed. In particular, we just modify the standard MESFET 
setup by allowing for ferromagnetic regions. For a more detailed modeling, we refer e.g. to 

na. 

In the following, we describe the geometry of the device in the {x,y) plane (see Figure 
[^. The total length is L = 0.6/im and the height equals H = 0.2/im. The source and 
drain regions are highly doped with doping = 3 ■ lO^^rn”^. The doping in the channel 
region is Cq = lO^^rn”^. The length of the source and drain regions are £ = 0.1 /xm. The 
gate contacts are attached at the middle of the device with a length of Lq = 0.2 /xm. The 
electrical parameters are: diffusion coefficient D = 10“^m^s“^, relaxation time r = Ips, 
temperature T = 300 K, and relative permittivity of silicon Sr = 11.7. These parameters 
are similar to those used in [13] (there is a small difference in the relaxation time value). 
With these data, the scaled Debye length is Xd = (eoErksT)/{qC+L"^) ~ 1.6 ■ 10““^, where 
Eq ~ 8.9 ■ 10“^^Fm“^ is the permittivity of the vaccum, q k. 1.6 ■ 10“^° C is the electron 
charge, and ks ~ 1.4 ■ 10“^^ m^kgs“^K“^ is the Boltzmann constant. 
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Figure 1. Geometry of a MESFET with ferromagnetic (F) source and drain 
regions and nonmagnetic (N) channel region. 


The gate contact is considered as a Schottky contact with barrier potential Vs = 0.8 V. 
The total voltage between source and gate is Vd + Vs, where Vg is the voltage applied at 
the gate. The density boundary value at the gate contact is calclulated according to PH 
Formula (5.1-19)], and the potential of the closed state is taken from [T3]. This gives 

• at the source: Rq = C+, n = 0, potential: 0 V, 

• at the drain: Uq = C+, n = 0, potential: Vd, 

• at the gate: 

open state: Uq = 3.9 ■ 10^^ m“^, n = 0, potential: Vs, 
closed state: Uq = 3.2 ■ 10® m“^, n = 0, potential: -|- 1.2 V, 

• for the other segments: homogeneous Neumann boundary conditions. 

The magnetic held is caused by the local orientation of the electron spin in the crystal and 
is predetermined by the ferromagnetic properties of the material. We consider a constant 
magnetic held, oriented along the ^-axis (perpendicular to the device). The electron spin 
may be also changed under the inhuence of the spin current, but we do not consider this 
ehect here. In our model, m corresponds to the direction of the local magnetic held, and 
the parameter 7 describes the intensity of the spin precession around this held. We choose 
m = 0 in the channel region and 


J (0, 0,1) for X < L/3 or a; > 2L/3, 
\ (0, 0, 0) for L/3<x < 2L/3. 


The value for 7 is taken from (TH], i.e. 7 = h/r, with h being the reduced Planck constant. 
The spin polarization is nonzero only in the highly doped source and drain regions, and 
we take p = 0.9. 

For the numerical discretization, we have chosen an admissible triangular mesh. Equa¬ 
tions 0 -(§ are approximated by scheme ([20|)-([25|), with the corresponding boundary 
conditions. The nonlinear system is solved at each time step by Newton’s method. The 
time step size is At = 0.05. The computations are continued until a steady state is reached 
or, more precisely, until the diherence of the solutions at two consecutive time steps in the 

norm falls below a threshold (typically, 10 “®). 
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Figure [^illustrates the scaled steady-state charge density no and the spin density na (note 
that ni = n 2 = 0) in the open state. The densities are scaled by the doping concentration 
C+, the spatial variable by the device length L. Compared to the closed state in Figure 
the charge density is rather large in the channel region, which can be also observed in 
standard MESFET devices. The charge current density in the closed state is by the factor 
10®... 10® smaller than in the open state. The spin density is (almost) zero in the closed 
state. 



0 0 0 0 


Figure 2. Scaled stationary charge density (left) and spin density (right) 
in an open-state MESFET with Vd = —2 V and Vg = 0 V. 

Current-voltage characteristics for MESFETs with and without ferromagnetic regions 
are shown in Figure [^ We observe that in the open state (nonpositive gate potentials), the 
current densities in the ferromagnetic MESFET are slightly larger than in the standard 
device, which allows for an improved device performance. When the transistor is closed 
(Vg = 1.2 V), the current densities are (almost) zero for both transistor types. 

In Figure [^ (left), we present the transient behavior of the charge density when switching 
from the open to the closed state (Vo = —2 V). The current values stabilize after about 1 ps. 
This justihes to define the numerical solution after 12ps as the “steady-state solution”. 
We compare these values with those computed from a standard MESFET; see Figure 
(right). The stabilization in the ferromagnetic case is slightly faster which allows for faster 
devices. 

Finally, we illustrate the free energy decay in Figure In this experiment, we have set 
Vd = 0 (source-drain voltage) and Vg = 0 (source-gate voltage). It turns out that the 
free energy decays with an exponential rate. For times larger than about 18 ps, the steady 
state is almost attained, and we observe numerical oscillations caused by the finite machine 
precision. 
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Figure 3. Scaled stationary charge density in a closed-state MESFET with 
Vd = —2V and Vq = 1.2V. 
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Figure 4. Current-voltage characteristics for the ferromagnetic (FM) and 
standard (NM) MESFET for various gate voltages Vq- For convenience, the 
source-drain voltages are given by their absolute values. 
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